;+
; NAME:
;       FILL_IN_RB3_GRS
;
; PURPOSE:
;       Fill in the rb3 structure with relevant IRDC model fit and GRS
;       information needed for later analyses.
;
; CATEGORY:
;       IRDC Distance Analysis
;
; CALLING SEQUENCE:
;       FILL_IN_RB3_GRS
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CASENO -- Case number to use.  [Default: 4].
;       ANGSEP -- Angular separation cutoff for GRS-association
;                 [arcmin].  [Default: 4]
;       SUFF   -- File suffix to use (i.e. '_q8').  [Default: '']
;
; KEYWORD PARAMETERS:
;       CLEM   -- Use the Clemens (1985) rotation curve for
;                 calculation of model fit distances from N/F fit
;                 placements.
;
; OUTPUTS:
;       NONE
;
; OPTIONAL OUTPUTS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  07/04/12, TPEB -- Code from check_grs_match.pro,
;                                   consolidated into a routine with
;                                   optional inputs.
;
;-

PRO FILL_IN_RB3_GRS, caseno, ANGSEP=angsep, CLEM=clem, SUFF=suff
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Parse inputs
  IF n_elements(caseno) EQ 0 THEN caseno=2
  IF caseno LT 0 || caseno GT 3 THEN message,'Improper model case number!'
  cno = string(caseno,format="(I1)")
  rotc = KEYWORD_SET(clem) ? 2b : 0b
  angsep = (n_elements(angsep) EQ 0) ? 4. : float(angsep)
  IF n_elements(suff) EQ 0 THEN suff = ''
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
  
  ;; Load everything
  IF n_elements(rb3) EQ 0 THEN $
     restore,'./irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  restore,'./ancillary/grs_kda_data.sav',/ver
  restore,'./irdc_dist_model/data/case'+cno+'_irdc_TdFIT'+suff+'.sav',/ver
  n = n_elements(rb3)
  
  rb3.constrain = constrain
  
  ;; Create arrays used for GRS identification
  dsep = fltarr(n)
  vsep = fltarr(n)
  grsi = lonarr(n)
  
  ;; Loop through RB3 to check for GRS co-incidence
  message,'Checking for GRS co-incidence...',/inf
  FOR i=0L, n-1 DO BEGIN
     gcirc, 2, rb3[i].l, rb3[i].b, grs.glon, grs.glat, dis
     mind = min(dis,gind)
     vsep[i] = abs(rb3[i].vlsr - grs[gind].vlsr)
     
     ;; A DSEP =-1 means there is no coincidence for this object
     dsep[i] = (vsep[i] LE grs[gind].delv*1.0) ? mind/60. : -1
     IF dsep[i] GT angsep THEN dsep[i] = -1
     grsi[i] = gind
  ENDFOR
  
  ;; Clean certain rb3 elements
  rb3.grs = 0b
  rb3.fitdist = 0.
  
  ;; Load rb3
  rb3.grs_clemdist = grs[grsi].dis
  rb3.grs_v        = grs[grsi].vlsr
  rb3.grs_l        = grs[grsi].glon
  rb3.grs_b        = grs[grsi].glat
  
  ;; Derive GRS distance based on Reid (2009) rotation curve
  grs_dtan = !MW.R0 * cos(rb3.grs_l * !dtor) / cos(rb3.grs_b * !dtor)/1.d3
  grs_far = rb3.grs_clemdist GE grs_dtan
  print,'N GRS at far KDIST: ',long(total(grs_far))
  
  
  ;; Check for match
  rb3.grs = dsep GE 0.
  message,'There are '+string(total(rb3.grs),format="(I0)")+$
          ' GRS-associated sources.',/inf
  
  ;; Read in GRS distance correction file
  readcol,'./irdc_dist_model/bad_grs_distances.txt',bad_gd_cnum,new_gd,$
          format='I,F'
  bgdind = WHERE_ARRAY(bad_gd_cnum,rb3.cnum,nbgd)
  rb3[bgdind].grs_dist = new_gd

  ;;==================================================================
  ;; Now, calculate whether a source is a GRS-match (based on the model)!
  
  ;; Go through each source and do the KDIST
  FOR i=0L, n-1 DO BEGIN
     
     IF (i+1) MOD 100 EQ 0 THEN $
        print,'Processing object '+string(i+1,format="(I3)")+$
              ' of '+string(n,format="(I0)")

     ;; Calculate grs_dist based on GRS KDA resolution
     rb3[i].grs_dist=(grs_far[i]) ? $
                     KDIST(rb3[i].grs_l,rb3[i].grs_b,rb3[i].grs_v,/FAR)/1.d3:$
                     KDIST(rb3[i].grs_l,rb3[i].grs_b,rb3[i].grs_v,/NEAR)/1.d3
     
     ;; Calculate the FITDIST based on the model N/F placement
     CASE rotc OF
        0: rb3[i].fitdist=(nearfar[i]) ? $
                          KDIST(rb3[i].l,rb3[i].b,rb3[i].vlsr,/FAR)/1.d3:$
                          KDIST(rb3[i].l,rb3[i].b,rb3[i].vlsr,/NEAR)/1.d3
        
        2: rb3[i].fitdist=(nearfar[i]) ? $
                          KDIST_CLEM(rb3[i].l,rb3[i].b,rb3[i].vlsr,/FAR)/1.d3:$
                          KDIST_CLEM(rb3[i].l,rb3[i].b,rb3[i].vlsr,/NEAR)/1.d3
        
     ENDCASE
     
     ;; GRS-match if delta_d <= 1 kpc
     rb3[i].grs_match = (abs(rb3[i].grs_dist - rb3[i].fitdist) LE 1.) ? 1b : 0b
     
  ENDFOR
  
  ;; Save rb3 file again
  save,rb3,filename='./irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  
END
